


###### INITIALIZE

setwd("C:/Users/spadmin/PhD/Projects/EPREmergence/Stats/SAUS/SAUS-NNM v1.0")
library(sandwich)
library(lmtest)
library(AER)
source("clusterrobust.R")


###### FUNCTIONS

deleteOverlaps <- function(input.df) {
  # Get points without overlaps
  single.points <- unique(input.df$rpid[is.na(input.df$rpid2)])
  
  # Iterate through points with overlaps and randomly eliminate overlapping points (i.e., buffers)
  input.op.df <- input.df[!is.na(input.df$rpid2), names(input.df) %in% c("rpid", "rpid2")]
  # Shuffle data frame
  input.op.df <- input.op.df[sample(1:nrow(input.op.df)),]
  # Iterate through overlapping points and eliminate overlaps
  op.cand.points <- unique(input.op.df$rpid)
  op.keep.points <- op.cand.points
  for (i in 1:length(op.cand.points)){
    thispoint <- op.cand.points[i]
    if (thispoint %in% op.keep.points){
      oppoints <- input.op.df$rpid2[input.op.df$rpid == thispoint]
      op.keep.points <- op.keep.points[!(op.keep.points %in% oppoints)]
    }
  }
  
  return(c(op.keep.points, single.points))
}


###### GET AND MODIFY BUFFER SAMPLE

# Get data from DB
rb.df <- dbGetTable(con, "hunzikp", "randombuffersdata")

# Data manipulation
rb.df$lnpop <- log(rb.df$pop1990 + 1)
rb.df$lnarea <- log(rb.df$area_sqkm)
rb.df$lncap <- log(rb.df$capdist_km + 1)
rb.df$lnborder <- log(rb.df$border_km + 1)
rb.df$cowid <- as.factor(rb.df$countries_cowid)
rb.df$groupcount09[is.na(rb.df$groupcount09)] <- 0
rb.df$groupcount65[is.na(rb.df$groupcount65)] <- 0
rb.df$gpdiff <- rb.df$groupcount09 - rb.df$groupcount65
rb.df$mgs09 <- ifelse(is.na(rb.df$min_groupsize09), 1, rb.df$min_groupsize09)
rb.df$max_noauton09 <- ifelse(rb.df$max_discrim09 == 1 | rb.df$max_powerless09==1, 1, 0)
rb.df$maxnoauton09[is.na(rb.df$max_noauton09)] <- 0
rb.df$max_excl09[is.na(rb.df$max_excl09)] <- 0
rb.df$max_excl[is.na(rb.df$max_excl)] <- 0
rb.df$max_powerless[is.na(rb.df$max_powerless)] <- 0
rb.df$max_discrim[is.na(rb.df$max_discrim)] <- 0
rb.df$mno <- ifelse(rb.df$max_powerless + rb.df$max_discrim > 0, 1, 0)
rb.df$lec <- log(rb.df$ethnologue_count+1)

# Subsetting
rb.df <- subset(rb.df, (ssafrica == 1 | asia == 1))  # Regional subsetting
rb.df <- subset(rb.df, cntr_petropoint == 1)  # Only countries with petroleum fields
rb.df <- rb.df[!is.na(rb.df$lnpop),]
rb.df <- rb.df[!is.na(rb.df$lnborder),]




###### RANDOM THINNING WITHOUT MATCHING

# Get table with overlaps for petropoints
op.df <- dbRunScript(con, "getBufferOverlaps_v0.1.sql", return=TRUE)

# Get IDs for randomly thinned sample
no.points <- deleteOverlaps(op.df)

# Get randomly thinned data frame
th.df <- rb.df[rb.df$rpid %in% no.points,]



###### ANALYSIS

## VALIDITY OF INSTRUMENT
pet.fit <- glm(petropoint ~ sb, data=th.df, family="binomial")
summary(pet.fit)
pet.fit <- glm(petropoint ~ sb + cowid + lnarea + lnpop + lncap + lnborder + elevsd + log(ethnologue_count+1), data=th.df, family="binomial")
summary(pet.fit)


## IVREG EXCLUSION
iv.fit <- ivreg(max_excl09 ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec | sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

iv.fit <- ivreg(max_excl09 ~ cowid + petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec | cowid + sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

iv.fit <- ivreg(mean_excl ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec | sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

iv.fit <- ivreg(mean_excl ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec + cowid | cowid + sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

# CLUSTER ROBUST S.E.s
# Asymptotic calculation
coeftest.cluster(th.df, iv.fit, cluster1="countries_cowid") 

# Cluster Bootstrap
iv.boot <- clBoot(iv.fit, th.df, "countries_cowid", reps=250, subst.coef=7, verbose=FALSE)







## IVREG NOAUTON
iv.fit <- ivreg(mean_noauton ~ petropoint | sb , data=th.df)
summary(iv.fit)

iv.fit <- ivreg(mean_noauton ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec | sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

iv.fit <- ivreg(mean_noauton ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec + cowid | cowid + sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

# CLUSTER ROBUST S.E.s
# Asymptotic calculation
coeftest.cluster(th.df, iv.fit, cluster1="countries_cowid") 

# Cluster Bootstrap
iv.boot <- clBoot(iv.fit, th.df, "countries_cowid", reps=250, subst.coef=7, verbose=FALSE)
iv.boot <- clBoot(iv.fit, th.df, "countries_cowid", reps=250, verbose=FALSE)





## IVREG GROUP COUNT
iv.fit <- ivreg(groupcount09 ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec | sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

iv.fit <- ivreg(groupcount09 ~ petropoint + lnpop + lnarea + lncap + lnborder + elevsd + lec + cowid | cowid + sb + lnpop + lnarea + lncap + lnborder + elevsd + lec, data=th.df)
summary(iv.fit)

# CLUSTER ROBUST S.E.s
# Asymptotic calculation
coeftest.cluster(th.df, iv.fit, cluster1="countries_cowid") 

# Cluster Bootstrap
iv.boot <- clBoot(iv.fit, th.df, "countries_cowid", reps=250, subst.coef=7, verbose=FALSE)